---
title: "List experiment Palestine"
author: "Dai Yamao"
date: "`r format(Sys.time(), '%Y-%m-%d')`"
output: html_document
---

```{r setup, include=FALSE}
rm(list=ls())
knitr::opts_chunk$set(echo = TRUE, fig.width = 10, fig.height = 5)
library(haven)
require(memisc)
require(stringi)
library(stargazer)
library(coefplot)
library(interplot)
library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(stargazer)
library(stringr)
library(tidyverse)
library(ggplot2)
library(psych)
#library(margins)
library(dplyr)
library(modelsummary)
library(makedummies)
library(ggeffects)
library(estimatr)
library(ggpubr)
library(BalanceR)

theme_set(theme_sjplot())

# PL
dat_pl <- read_csv("data/Palestine.csv")
dat_pl <- dat_pl[-1,]

dat_pl <- dat_pl |>
  mutate(treatment = `Q12-t_1`) |>
  mutate(control = `Q12-c_1`)|>
  mutate(List_experiment_Hamas = case_when(control >= 0 ~  "Control",
                                           treatment >= 0 ~ "Treatment")) |>
  mutate(Treatment = ifelse(List_experiment_Hamas == "Treatment", "1", "0"))

dat_pl$`Q12-c_1`[is.na(dat_pl$`Q12-c_1`)] <- 0
dat_pl$`Q12-t_1`[is.na(dat_pl$`Q12-t_1`)] <- 0

dat_pl <- dat_pl |>
  mutate(Hamas_experiment = `Q12-c_1` + `Q12-t_1`)

dat_pl <- dat_pl |>
  mutate(Gender = D1,
         Age = D2,
         Education = D3,
         Income = D5)

# LB
dat_lb <- read_csv("data/Lebanon.csv")
dat_lb <- dat_lb[-1,]

dat_lb <- dat_lb |>
  mutate(treatment = `Q12-t_1`) |>
  mutate(control = `Q12-c_1`)|>
  mutate(List_experiment_Hezbollah = case_when(control >= 0 ~  "Control",
                                               treatment >= 0 ~ "Treatment")) |>
  mutate(Treatment = ifelse(List_experiment_Hezbollah == "Treatment", "1", "0"))

dat_lb$`Q12-c_1`[is.na(dat_lb$`Q12-c_1`)] <- 0
dat_lb$`Q12-t_1`[is.na(dat_lb$`Q12-t_1`)] <- 0

dat_lb <- dat_lb |>
  mutate(Hezbollah_experiment = `Q12-c_1` + `Q12-t_1`)

dat_lb <- dat_lb |>
  mutate(Gender = D1,
         Age = D2,
         Education = D3,
         Income = D5)

# IQ
dat_iq <- read_csv("data/Iraq.csv")
dat_iq <- dat_iq[-1,]

dat_iq <- dat_iq |>
  mutate(treatment = `Q12-t_1`) |>
  mutate(control = `Q12-c_1`)|>
  mutate(List_experiment_PMU = case_when(control >= 0 ~  "Control",
                                         treatment >= 0 ~ "Treatment")) |>
  mutate(Treatment = ifelse(List_experiment_PMU == "Treatment", "1", "0"))

dat_iq$`Q12-c_1`[is.na(dat_iq$`Q12-c_1`)] <- 0
dat_iq$`Q12-t_1`[is.na(dat_iq$`Q12-t_1`)] <- 0

dat_iq <- dat_iq |>
  mutate(PMU_experiment = `Q12-c_1` + `Q12-t_1`)

dat_iq <- dat_iq |>
  mutate(Gender = D1,
         Age = D2,
         Education = D3,
         Income = D5)

# YM
dat_ym <- read_csv("data/Yemen.csv")
dat_ym <- dat_ym[-1,]

dat_ym <- dat_ym |>
  mutate(treatment = `Q12-t_1`) |>
  mutate(control = `Q12-c_1`)|>
  mutate(List_experiment_Houthi = case_when(control >= 0 ~  "Control",
                                            treatment >= 0 ~ "Treatment")) |>
  mutate(Treatment = ifelse(List_experiment_Houthi == "Treatment", "1", "0"))

dat_ym$`Q12-c_1`[is.na(dat_ym$`Q12-c_1`)] <- 0
dat_ym$`Q12-t_1`[is.na(dat_ym$`Q12-t_1`)] <- 0

dat_ym <- dat_ym |>
  mutate(Houthi_experiment = `Q12-c_1` + `Q12-t_1`)

dat_ym <- dat_ym |>
  mutate(Gender = D1,
         Age = D2,
         Education = D3,
         Income = D5)
```

# balance check
```{r}
# PL
balance_pl <- BalanceR(data = dat_pl, 
                    group = List_experiment_Hamas,
                    cov = c(Gender, Age, Education, Income))

balance_PL <- plot(balance_pl, point.size = 5, text.size = 18) +
  labs(title = "Palestine")
balance_PL

# LB
balance_lb <- BalanceR(data = dat_lb, 
                    group = List_experiment_Hezbollah,
                    cov = c(Gender, Age, Education, Income))

balance_LB <- plot(balance_lb, point.size = 5, text.size = 18) +
  labs(title = "Lebanon")
balance_LB

# IQ
balance_iq <- BalanceR(data = dat_iq, 
                    group = List_experiment_PMU,
                    cov = c(Gender, Age, Education, Income))
balance_IQ <- plot(balance_iq, point.size = 5, text.size = 18) +
  labs(title = "Iraq: List experiment")
balance_IQ

# YM
balance_ym <- BalanceR(data = dat_ym, 
                    group = List_experiment_Houthi,
                    cov = c(Gender, Age, Education, Income))
balance_YM <- plot(balance_ym, point.size = 5, text.size = 18) +
  labs(title = "Yemen: List experiment")
balance_YM

# plot
p10 <- ggarrange(balance_PL,
                 balance_LB + theme(axis.title.y = element_blank()),
                 balance_IQ,
                 balance_YM + theme(axis.title.y = element_blank()),
                 nrow = 2, ncol = 2, align = "hv")
p10
ggsave(filename = "result/balance_list_1.png",
       plot = p10,
       width = 15,
       height = 8,
       units = "in", 
       dpi = 600)
```

# List results
```{r}
# PL
list_pl <- dat_pl |> 
  summarise(control_mean = mean(control, na.rm = TRUE), control_SE = sd(dat_pl$control, na.rm = TRUE) / sqrt(length(dat_pl$control)),
            treatment_mean = mean(treatment, na.rm = TRUE), treatment_SE = sd(dat_pl$treatment, na.rm = TRUE) / sqrt(length(dat_pl$treatment)))

df_pl <- data.frame(
  Experiment = c("Control", "Treatment"),
  Mean = c(list_pl$control_mean, list_pl$treatment_mean),
  SE = c(list_pl$control_SE, list_pl$treatment_SE)
)

pl <- ggplot(df_pl, aes(x = Experiment, y = Mean, group = Experiment)) +
  geom_point(stat = "identity", color = "black") +
  geom_errorbar(aes(ymax = Mean+SE, ymin = Mean-SE), width = 0.2, position = position_dodge(.9)) +
  labs(title = "List expriment of Hamas") +
  theme_bw()
pl

# LB
list_lb <- dat_lb |> 
  summarise(control_mean = mean(control, na.rm = TRUE), control_SE = sd(dat_lb$control, na.rm = TRUE) / sqrt(length(dat_lb$control)),
            treatment_mean = mean(treatment, na.rm = TRUE), treatment_SE = sd(dat_lb$treatment, na.rm = TRUE) / sqrt(length(dat_lb$treatment)))

df_lb <- data.frame(
  Experiment = c("Control", "Treatment"),
  Mean = c(list_lb$control_mean, list_lb$treatment_mean),
  SE = c(list_lb$control_SE, list_lb$treatment_SE)
)

lb <- ggplot(df_lb, aes(x = Experiment, y = Mean, group = Experiment)) +
  geom_point(stat = "identity", color = "black") +
  geom_errorbar(aes(ymax = Mean+SE, ymin = Mean-SE), width = 0.2, position = position_dodge(.9)) +
  labs(title = "List expriment of Hezbollah") +
  theme_bw()
lb

# IQ
list_iq <- dat_iq |> 
  summarise(control_mean = mean(control, na.rm = TRUE), control_SE = sd(dat_iq$control, na.rm = TRUE) / sqrt(length(dat_iq$control)),
            treatment_mean = mean(treatment, na.rm = TRUE), treatment_SE = sd(dat_iq$treatment, na.rm = TRUE) / sqrt(length(dat_iq$treatment)))

df_iq <- data.frame(
  Experiment = c("Control", "Treatment"),
  Mean = c(list_iq$control_mean, list_iq$treatment_mean),
  SE = c(list_iq$control_SE, list_iq$treatment_SE)
)

iq <- ggplot(df_iq, aes(x = Experiment, y = Mean, group = Experiment)) +
  geom_point(stat = "identity", color = "black") +
  geom_errorbar(aes(ymax = Mean+SE, ymin = Mean-SE), width = 0.2, position = position_dodge(.9)) +
  labs(title = "List expriment of PMU") +
  theme_bw()
iq

# YM
list_ym <- dat_ym |> 
  summarise(control_mean = mean(control, na.rm = TRUE), control_SE = sd(dat_ym$control, na.rm = TRUE) / sqrt(length(dat_ym$control)),
            treatment_mean = mean(treatment, na.rm = TRUE), treatment_SE = sd(dat_ym$treatment, na.rm = TRUE) / sqrt(length(dat_ym$treatment)))

df_ym <- data.frame(
  Experiment = c("Control", "Treatment"),
  Mean = c(list_ym$control_mean, list_ym$treatment_mean),
  SE = c(list_ym$control_SE, list_ym$treatment_SE)
)

ym <- ggplot(df_iq, aes(x = Experiment, y = Mean, group = Experiment)) +
  geom_point(stat = "identity", color = "black") +
  geom_errorbar(aes(ymax = Mean+SE, ymin = Mean-SE), width = 0.2, position = position_dodge(.9)) +
  labs(title = "List expriment of Houthis") +
  theme_bw()
ym

# plot
p11 <- ggarrange(pl + theme(axis.title.x = element_blank()),
                 lb + theme(axis.title.x = element_blank()) + theme(axis.title.y = element_blank()),
                 iq,
                 ym + theme(axis.title.y = element_blank()),
                 nrow = 2, ncol = 2, align = "hv")
p11
ggsave(filename = "result/list1_result.png",
       plot = p11,
       width = 9,
       height = 7,
       units = "in", 
       dpi = 600)
```